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ABSTRACT 

We carry out local three dimensional (3D) hydrodynamic simulations of planet-disk interaction in 
stratified disks with varied thermodynamic properties. We find that whenever the Brunt- Vaisala 
frequency (N) in the disk is nonzero, the planet exerts a strong torque on the disk in the vicinity of 
the planet, with a reduction in the traditional "torque cutoff". In particular, this is true for adiabatic 
perturbations in disks with isothermal density structure, as should be typical for centrally irradiated 
protoplanetary disks. We identify this torque with buoyancy waves, which are excited (when N is 
non-zero) close to the planet, within one disk scale height from its orbit. These waves give rise to 
density perturbations with a characteristic 3D spatial pattern which is in close agreement with the 
linear dispersion relation for buoyancy waves. The torque due to these waves can amount to as much 
as several tens of per cent of the total planetary torque, which is not expected based on analytical 
calculations limited to axisymmetric or low-m modes. Buoyancy waves should be ubiquitous around 
planets in the inner, dense regions of protoplanetary disks, where they might possibly affect planet 
migration. 

Subject headings: hydrodynamics, waves, stars: formation, stars: pre-main sequence, planet-disk in- 
teractions 



1. INTRODUCTION 

The gravitational potential of a planet surrounded 
by a protoplanetary disk is known to give rise to non- 
axisymmetric density waves, which propagate away from 
the planet carrying angular momentum. Gravitational 
coupling with these density perturbations exerts a torque 
on the planet, which might lead to its migration if the 
one-sided torques produced by the inner and outer disks 
do not cancel. Pioneering linear calculations by Goldre- 
ich & Tremaine (1980) have shown that, neglecting this 
possible small asymmetry, the one-sided torque acting on 
a planet moving on a circular orbit in a razor-thin (2D) 
disk is 



T = q 2U (GM p 



2 ^()RpQp 



? 2D « 0.93, 



(1) 



where M P , R p and £l p are the planetary mass, semi- 
major axis and angular frequency respectively, c s = 
(dP/dpls) 1 ^ 2 is the adiabatic sound speed, and £o is the 
unperturbed disk surface density. Most of this torque 
is excited at Lindblad resonances corresponding to the 
m ~ R p /h 3> 1 azimuthal harmonic of the planetary po- 
tential and located at a distance ~ h — c s /tt p away from 
the planetary orbit. Closer to the planet, at separations 
< h the excitation torque density decreases exponentially 
— the so-called "torque cut off" phenomenon. 

In three-dimensional (3D) stratified disks a variety of 
waves can be excited (Lubow & Pringle 1993, Korycan- 
sky & Pringle 1995, Lubow & Ogilvie 1998), including 
buoyancy waves whenever the Brunt- Vaisala frequency 
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is non-zero. Global analytical calculation of the wave 
properties for all modes are available only for one spe- 
cial case, in which both the vertical disk structure and 
the equation of state (EOS) are isothermal, and N = 0. 
In this case, the result (1) still holds but with 52D low- 
ered to <73D ~ 0.25(?2D = 0.37 as a result of averaging 
the planetary potential over the vertical disk structure 
(Takeuchi & Miyama 1998). For other combinations of 
the vertical disk structure and EOS, in particular those 
characterized by non-zero TV, analytical results are lim- 
ited only to either axisymmetric or low-m modes excited 
far from the planet. In that case Lubow & Ogilvie (1998) 
found that although buoyancy waves are responsible for 
a fraction of the total torque excited by the planet, the 
total torque from all modes is the same as that given by 
traditional Lindblad torque formulae. 

In this paper we carry out 3D simulations of stratified 
disks characterized by non-zero N with embedded low- 
mass planets to explore global excitation of buoyancy 
waves and their role in angular momentum exchange with 
the planet. 

2. NUMERICAL METHOD 

The numerical tool we used is Athena (Stone et al. 
2008), a grid-based code with a higher-order Godunov 
scheme, piecewise parabolic method (PPM) for spatial 
reconstruction, and the corner transport upwind (CTU) 
method for multidimensional integration. We use the 
stratified shearing-box set-up implemented in Athena 
(Stone & Gardiner 2010). Since the potential vortic- 
ity is zero in the shearing-box, the traditional corotation 
torque is zero (Goldreich & Tremaine, 1979). 

A planet is placed at the origin of a Cartesian box in 
x,y,z coordinates. Its smoothed potential is approxi- 
mated as 



$ p (d) = -GM P 
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where r s is the smoothing length and d — \/ x 2 + y 2 + z 2 . 
This potential converges to the point mass potential as 
(r s /d) 4 for d >• r s . In most runs r s = 0.125/i, which is 
resolved by 4 grids cells with our normal resolution of 
32/h. 

To reduce the computation time we take advantage of 
the symmetry of the problem and only simulate the disk 
for x = [-X,0], y = [-Y, Y], and z = [0,Z], see Figure 
1. We take X = 5h, Y = 20ft. so that the planetary wake 
located at y <~ — 3x 2 /4/i fits well inside the box. We vary 
Z in different runs (see Table 1) in such a way that the 
density at z — Z is ^5 orders of magnitude lower than 
the midplane density. 

We use different boundary conditions (BCs) at differ- 
ent box faces: outflow BC at the y = Y face, and "fixed 
state" BC (meaning that all physical variables in the 
ghost zones are fixed at their unperturbed Keplerian val- 
ues) at x = —X and y = —Y faces (for one model we also 
use periodic BC in y-direction, see Table 1). At the x = 
face we employ a "symmetric" BC: the x = face has 
been divided into two portions — y > and y < 0, and 
the variables in the last 4 active zones of one portion are 
copied into the ghost zones in the other portion symmet- 
ric with respect to x = y = 0; velocities (or momenta) in 
x and y directions are copied with the opposite sign. We 
use a reflecting BC at z — since both the planetary po- 
tential and the disk structure are symmetric with respect 
to z. At z = Z we use an outflow BC but the densities 
are extrapolated assuming hydrostatic equilibrium into 
the ghost zones, and to prevent inflow if v z < 0, then we 
set v z to 0. 

2.1. Physical setup 

In our calculations we use both an isothermal (7=1 
in p oc p 7 ), and an adiabatic EOS (7 = 5/3 or 7/5). 
The unit of length in our simulations is defined via the 
adiabatic scale height h = c s /Q p = ^/7~Cj so /f2, where 

c iso = ikT/p) 1 ! 2 and c s is the adiabatic sound speed 
calculated using the midplane temperature. Obviously, 
c s = Ci SO for an isothermal EOS (7 = 1). This choice 
keeps sound waves traveling the same distance in a given 
amount of time in different models so that the wake po- 
sition is always the same. 

The vertical structure of a protoplanetary disk is typ- 
ically determined by irradiation from the central star so 
that is vertically isothermal 
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where poo is the midplane density. The value of Ci SO 
is determined by the central irradiation flux (Kcnyon & 
Hartmann 1987; Chiang & Goldreich 1997). 

We also consider thermally stratified disks in which 
the disk midplane is hotter than its atmosphere, which is 
expected if viscous heating is important. Including the 
fact that such disks will eventually become isothermal 
at hight z, the vertical structure of the disk is better 
represented as 



Po(z) = c 2 s , s po{z) 



Po{z) 



r-i 



+ 1 



(5) 



speed and density at the transition where the disk be- 
comes isothermal (Lin et al. 1990; Bate ct al. 2002). How- 
ever, unlike Bate et al. (2002), we do not necessarily set T 
in Eq. (5) the same as 7 in the EOS. We derive po{z) by 
solving the equation of vertical hydrostatic equilibrium 
with Eq. (5) using the Newton-Raphson method. When 
r — > 1 the isothermal disk structure is recovered. 

We run a set of five models up to 10 orbits with dif- 
ferent disk structure and EOS, summarized below and in 
Table 1. 

(a) Isothermal disk structure (4) with isothermal EOS 
(model II). 

(b) Disk structure (5) with adiabatic EOS and T = 
7 = 5/3 (model AA). 

(c) Isothermal disk structure with adiabatic EOS hav- 
ing 7 = 5/3 (model IA). This setup is typical for a cen- 
trally irradiated disk in which density perturbations are 
adiabatic. 

(d) Analogous to (4) but with 7 = 7/5 (model IA2), 
designed to illustrate the effect of adiabatic index on the 
buoyancy wave coupling. 

(e) Disk structure (5) with r = 5/3 and adiabatic EOS 
with 7 = 7/5 (model PA), to illustrate viscously heated 
accretion disks. 

In the first two models N — so that buoyancy waves 
are absent, while the last three models feature non-zero 
N and support buoyancy waves. All models were run at 
a resolution of 32/h, but we also run model (3) at higher 
resolution (64//i) to verify numerical convergence. 

The planet mass is M p — 0.0058M t ^, where M t h is the 
thermal mass 



M th = 



(6) 



where c s . s and p s ( set as 0.01 p c ) represent the sound 



Thus, disk-planet coupling is well inside the linear 
regime. Density wave generated by such a planet in a 
2D disk would shock at \x\ = 6h (Goodman & Rafikov 
2001), which is outside our box. 

3. RESULTS 
3.1. Without Buoyancy Waves 

First we show the results for models with no buoy- 
ancy waves (II and AA) . The spatial distributions of the 
excitation torque density dT/dx (the amount of torque 
excited per unit radial distance) for these two models 
are shown in Fig. 2(a) (b) as green and orange curves. 
One can see that dT/dx rapidly goes to zero for \x\ < h, 
clearly exhibiting the "torque cut-off" phenomenon (Gol- 
dreich & Tremaine 1980) in both models. 

For model II, our simulation agrees with the semi- 
analytical calculation by Takeuchi & Miyama (1998) re- 
markably well, and we confirm that Eq. (1) holds in this 
case with q3r>(H) = 0.25q2D- Although wave structures 
are different in model AA (Lubow & Pringle 1993). the 
distribution of dT/dx for this model is quite similar in 
shape to that of model II, but the amplitude is slightly 
different. 

3.2. With Buoyancy waves 

We next look at the behavior of dT/dx in models with 
non-zero N. The most important feature of these models, 
evident in Figure 2, is the weak torque cutoff near the 
planet — all show a significant torque contribution near 
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x = 0, including model IA2, which provides the best 
description for an irradiated protoplanetary disk. Since 
higher-m modes are excited at Lindblad resonances closer 
to the planet, the strong torque close to the planet also 
suggests the traditional torque cutoff in Fourier space 
(Ward 1997) is weaker when buoyancy waves are present. 

The one-sided torque T (Fig. 2(c) (d)) in models 
IA and IA2 is characterized by q^(IA) = 0.48 and 
g 3 D(J^42) = 0.41 respectively, which should be compared 
with qsD(II) = 0.34 for model II having the same verti- 
cal structure as IA and IA2 but different EOS. 

In order to identify the origin of this excess torque near 
the planet, we integrate the volume density of the torque 
along y direction 

d2T f°° A A < \ d< t> (7\ 

——= dydp(x,y,z) — . (7) 
dxdz J_ oc dy 

This quantity is shown in Fig. 3(a) for models II and IA. 
We see that the excess torque in model IA comes from 
the region z <~ 0.4 h and \x\ < h. In Fig. 3(d) we plot the 
density contours in the xy plane at z = OAh. In clear 
contrast with model II (Fig. 3(c)), the density distri- 
bution in model IA exhibits density fluctuations/ridges 
close to the planet which are different from the usual 
wake structure. They extend to y > and look like rays 
emanating from the origin. It is these density fluctua- 
tions that contribute to the excess torque. 

These fluctuations also have vertical structure, as 
shown in Fig. 4 where we plot Sp = p — po and v z in the 
yz plane at x = 0.5ft for both models II and IA. Again, 
significant density fluctuations close to the planet exist 
only in case IA. 

Since the only difference between case IA and case II 
is that the case IA has non-zero Brunt- Vaisala frequency 
N, it is natural to relate these density fluctuations to 
buoyancy waves. To confirm this hypothesis, we note 
that a mode with wavenumber k y has the buoyancy fre- 
quency N whenever the condition 

2Ak y x = N (8) 

is fulfilled. For comparison, the Lindblad resonance con- 
dition is 2Ak y x = ±k. In the disk with isothermal struc- 
ture and adiabatic EOS (models IA and IA2), we have 

With 7 = 5/3, and X y — 2n/k y , Eqs. 8 and 9 can be 
combined to derive 

^ - 11.543-. (10) 
ft z 

If we assume the phase of buoyancy waves is at x=0, 
the geometric location of the constant phase 2nn (n is 
integer) is given by 

\ = 11.543-ra and n= 1,2,3,.... (11) 

ft z 

Curves corresponding to Eq. 11 are drawn in the lower 
right panels of Figs. 3 and 4. They agree with the pattern 
of the density perturbation and v z quite well, strongly 
suggesting that these fluctuations and the excess torque 
are related to buoyancy waves (their dissipation or exci- 
tation). 



The torque density and integrated torque for model 
IA2 (isothermal disk with adiabatic EOS and 7 = 1.4) 
plotted as the cyan curves in Fig. 2 also exhibit the ex- 
cess torque density due to buoyancy waves close to the 
planet. However, the integrated torque in this case is 
characterized by q$D {IA2) s=s 0.41, which is lower than 
qsD(IA), and the buoyancy wave coupling to the plane- 
tary potential is weaker for smaller 7 as seen from dT/dx 
curve. 

Blue curves in the right panels of Fig. 2 show dT/dx 
and T(x) for model PA. There is again an apparent 
torque excess near the planet, which is not surprising 
since N 7^ in this model allowing buoyancy waves to 
be excited. 

4. DISCUSSION 

The phenomenon of torque cut-off in 2D disks is due 
to the fact that high-fcj, modes (k y > ft -1 , excited close 
to the planet) of rotation-modified sound waves cou- 
ple poorly to the planetary potential (even though it is 
strongest there). However, the dispersion relation for 
buoyancy waves is very different from that of the waves 
in 2D disks. Based on Eq. (10), the pattern in our simu- 
lations should have k y w h~ l at x ~ 0.05ft- and z = O.l/i, 
resulting in strong coupling between buoyancy waves and 
the planetary potential. This may explain the lack of 
torque cutoff in models with non-zero N. A further in- 
vestigation will be presented in Zhu et al. (in prep.). 

Since the buoyancy torque is present very close to the 
planet, its strength could be sensitive to the potential 
softening length r s in the simulation. Thus, we re-ran 
model IA with a smaller smoothing length (r s — l/16/i) 
but obtained essentially the same result (purple curves in 
Fig. 2), suggesting that our numerical results are robust. 

4.1. Comparison with previous studies 

Previous analytical and semi-analytical studies of 
buoyancy waves (Lubow & Pringle 1993, Korycansky & 
Pringle 1995, Lubow & Ogilvie 1998) have been limited 
to axisymmctric waves or low-m modes excited far away 
from the planet, at \x\ » ft. These studies have demon- 
strated such modes to play rather insignificant role in 
planet-disk angular momentum exchange. 

The main result of this work is that coupling of the 
planetary potential to the locally- excited (\x\ < h) buoy- 
ancy waves channels a considerable amount of angular 
momentum into these modes. As a result, the buoyancy 
torque strongly contributes to the planetary torque, at 
the level of tens of per cent. This result is different from 
analytical expectations extrapolated from low-m buoy- 
ancy modes (Lubow & Ogilvie 1998). 

Most previous 3D simulations of planet-disk interac- 
tion explored isothermal disk structure with an isother- 
mal EOS (e.g. Bate et al. 2003, D'Angelo & Lubow 
2010). An adiabatic EOS has been explored in 
Paardekooper & Mellema (2006). However, because of 
the global geometry used in their work, the buoyancy 
torque may have been hard to distinguish from the coro- 
tation torque, see §4.3. 

4.2. Existence and Signatures of Buoyancy Waves 

The vertical structure of protoplanetary disks around 
T Tauri stars should be close to isothermal given the 
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dominant role of stellar irradiation in their thermal bal- 
ance. Existence of buoyancy waves in such disks de- 
pends on the thermodynamics of the fluid perturbations 
on timescales ~ iV _1 . Density perturbations induced by 
planets result in temperature fluctuations which tend to 
be erased by radiative diffusion. If the cooling time t c 
of such perturbations is longer than N^ 1 then modes 
are adiabatic with 7 > 1 (as in models IA and IA2) 
and buoyancy waves should be present. In the opposite 
case of t c N < 1 fluid perturbations behave as isothermal 
(7 w 1) and buoyancy waves are not excited, as in model 
II. 

We expect that the separation between the two regimes 
should occur between several to several tens of AU, de- 
pending on the disk properties, and that buoyancy waves 
should be efficiently excited by relatively close-in planets. 
Further out t c N < 1 and excitation of buoyancy waves 
should be suppressed. Note, that both t c and N are in 
general functions of z, sec equation (9). This makes cool- 
ing considerations dependent not only on the radius but 
also on the height in the disk. 

Whenever buoyancy waves are efficiently excited by 
planets they can reveal themselves via associated ver- 
tical motions (see Fig. 1) which should give rise to cor- 
rugations of the disk surface aligned with density ridges 
seen in Fig. 3. Ncar-IR imaging of stellar light scat- 
tered by such perturbations at very high resolution (on 
scales <~ h) can reveal even low-mass planets embedded 
in protoplanetary disks. 

4.3. Migration 

The rate at which a planet migrates is determined by 
the imbalance of one-sided torques exerted on it by the 
inner and outer parts of the disk, caused by the global ra- 
dial gradients of the disk density and temperature (Ward 
1986). The net torque can also be caused by the local 
asymmetries in the disk properties caused by the pres- 
ence of the planet itself, such as the temperature pertur- 
bations due to shadowing near the planet (Jang-Condell 
& Sasselov 2005), or the planet's own heat output caused 
by e.g. the accretion of planetesimals. Such asymmetries 
generated very close to the planet may not strongly affect 
the net standard Lindblad torque because of the torque 
cutoff at such separations. However, their effect on the 
net buoyancy torque excited primarily at small x may 
be disproportionately large and considerably affect the 
planetary migration. 

This statement is only true if the buoyancy waves ob- 
served in this work can propagate and deposit their angu- 
lar momentum far from the corotation region. Otherwise 
their angular momentum accumulates in the narrow an- 
nulus around the planetary orbit and the corresponding 
torque on the planet may be subject to saturation, like 
the standard corotation torque (Balmforth & Korycan- 
sky 2001; Masset 2001, 2002). This would eliminate the 
effect of the buoyancy waves on the planetary migration. 
We will investigate these possibilities in Zhu et al. (in 
prep). 

Direct measurement of the buoyancy torque effect on 
the migration speed would require global disk-planet cal- 
culations, in which corotation torque appears as well. It 
may then be difficult to separate the effects of the buoy- 
ancy and corotation torques in global simulations. How- 
ever, the strength of the corotation torque depends on the 



radial gradients of specific vortensity and entropy across 
the horseshoe region (Paardekooper & Mcllema 2006; 
Baruteau & Masset 2008; Paardekooper & Papaloizou 
2008). By properly setting the disk properties to nullify 
these gradients one can eliminate the corotation torque 
in global simulations, thus isolating the contribution due 
to the buoyancy torque. 

4.4. Wave dissipation and gap opening 

Gap opening by massive planets depends not only on 
the wave excitation but also on wave damping as a means 
of transferring angular momentum to the disk fluid (Lu- 
nine & Stevenson 1982; Rafikov 2002). The global den- 
sity wake excited by planet is thought to dissipate pri- 
marily via shock damping (Goodman & Rafikov 2001), 
but buoyancy waves are very distinct from rotation- 
modified sound waves and should dissipate differently 
(Lubow & Ogilvie 1998; Bate et al. 2002). Channeling 
of the wave action (Lubow & Ogilvie 1998) may con- 
siderably speed up their nonlinear evolution resulting in 
more efficient damping. Since we showed that buoyancy 
waves carry good fraction of the total angular momen- 
tum flux, understanding their damping mechanism and 
spatial pattern of dissipation may be important for clar- 
ifying the issue of gap opening by planets (Zhu et al. in 
prep). 

Our work suggests buoyancy waves can play an im- 
portant role in planet migration and gap opening which 
demands further studies. 

Authors are indebted to Jeremy Goodman, Steve 
Lubow, and Gordon Ogilvie for helpful comments and 
suggestions. This work was supported by NSF grant 
AST-0908269 and Princeton University. This research 
was supported in part by the NSF through TeraGrid 
resources provided by the Texas Advanced Computing 
Center and the National Institute for Computational Sci- 
ence under grant number TG-AST090106. 
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TABLE 1 
Models 



Case 


structure 


EOS 


y-direc. 


Brunt- Vaisala 


iorque 


Z domain size 


Resolution 






7 


boundary 


N 


coefficient (q) a 






XxYxZ 


II 


isothermal 


1 


in/ outflow 





0.34 




[0, 5 h] 


160x1280x160 


AA 


poly. T = 5/3 


5/3 


in/outflow 





0.44 




[0, 2 h] 


160x1280x64 


AApcr 


poly. T = 5/3 


5/3 


periodic 





0.44 




[0,2 h] 


160x1280x64 


IA 


isothermal 


5/3 


in/outflow 


7^0 


0.48 


[0, 3.873 h] 


160x1280x160 


IApcr 


isothermal 


5/3 


periodic 


7^0 


0.46 




0, 3.873h] 


160x1280x160 


IA (64/h) 


isothermal 


5/3 


in/outflow 


7^0 


0.48 




0, 3.873h] 


320x2560x320 


IAss (64/h, r s =l/16 h) 


isothermal 


5/3 


in/outflow 


7^0 


0.48 




0, 3.873h] 


320x2560x320 


IA2 


isothermal 


7/5 


in/outflow 


7^0 


0.41 




0, 3.873h] 


160x1280x160 


PA 


poly. T = 7/5 


5/3 


in/outflow 


7^0 


0.45 




[0, 2.5h] 


160x1280x80 



q as in T = q(GM p ) 2 Y, R v U/ c^ adi 
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Fig. 1. — The iso-density contour of p = 0.8poo for an IA disk model. Planetary mass is increased to M p = 0.3Mth to visually enhance 
the effect of buoyancy waves. Buoyancy waves cause ray-like density disturbances close to x = (right side of the box), which are also 
shown in Fig. 3. Two streamlines passing through (-0.5, 0, l)(blue) and (-3, 0, l)(red) demonstrate vertical oscillation due to buoyancy 
waves and the small velocity perturbation of sound waves. 
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Fig. 2. — Torque densities (dT/dx, upper panels) and integrated torques (lower panels) for disks with isothermal disk structure (left 
panels) and polytropic disk structure (right panels) with different EOS. Models If (green curves) and AA (orange curves) not supporting 
buoyancy waves exhibit a clear torque cut-off at \x\ < h. All other models support buoyancy waves and show significant buoyancy torque 
near the planet. 
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Fig. 3. — Upper panels: (volume) density of the planet-induced torque integrated over {/-coordinate for models II (left panel) and I A 
(right panel). Strong excess torque due to buoyancy waves close to the planet is clearly visible at z ~ h in model IA, which features non-zero 
Brunt- Vaisla frequency. Lower panels: density contours in the xy plane at z where po equal to 0.8755 poo (indicated by the dashed line 
at z ~ 0.5h in (a) and z ~ OAh in (b), h are different in these two cases). Geometric locations corresponding to the resonance condition 
for buoyancy waves (11) are shown by dotted lines in model IA, and agree quite well with the pattern of density fluctuations derived from 
simulations. 
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Fig. 4. — Upper panels: density fluctuations in the yz plane at \x\ = 0.5h for models II (left panel) and IA (right panel). Lower panels: 
vertical velocity v z at \x\ = 0.5h for models II and IA. The buoyancy resonance positions (11) are again plotted as dotted curves. The 
fluctuations of v z at z > 3h are due to the lack of exact hydrostatic equilibrium at the upper boundary. 



